---
title: "Climate Policy Beliefs, not Problem Visibility, Drive Rural Opposition to Wind Energy Projects: Replication"
author: Inhwan Ko
output: html_document
date: "2025-11-21"
---

```{r setup, include=FALSE}
#rm(list=ls())
pacman::p_load(MASS, readxl, tidyverse, sandwich, lmtest, mediation)
```

# 1. Prepare the data

## 1.1. Load the data

```{r}
data <- read_excel("rawdata.xlsx")
rid <- read.csv("rid.csv") %>% select(-X)
data <- left_join(rid, data, by="rid")

data <- data %>% filter(Progress==100) %>% # leave only completes
  select(-Progress, -RecordedDate, -'AC-Q1_DO', -'AC-Q2_DO', -'D-Q7_DO', -rid) 
  # unnecessary variables
colnames(data) <- c("startdate","enddate","dur","id",
                    "bornyear","zipcode","support","ac1","ac2","gender","race",
                    "reside","educ","income","occu","party","relig","urbanrural",
                    "q1","q2","q3","q4","q5","q6","q7","treat")
```

## 1.2. Filter out attention check failures 

```{r}
data$ac <- 0
data$ac[data$treat=="Treatment1" & # wind turbines in farmland
          data$ac1=="Wind turbines" & data$ac2=="Farmland with crops"] <- 1
data$ac[data$treat=="Treatment2" & # wind turbines in farmland
          data$ac1=="Wind turbines" & data$ac2=="Farmland with crops"] <- 1
data$ac[data$treat=="Treatment3" & # no turbines in farmland
          is.na(data$ac1) & data$ac2=="Farmland with crops"] <- 1
data$ac[data$treat=="Treatment4" & # wind turbines in industrial landscape
          data$ac1=="Wind turbines" & data$ac2=="Industrial buildings"] <- 1
data$ac[data$treat=="Treatment5" & # wind turbines in industrial landscape
          data$ac1=="Wind turbines" & data$ac2=="Industrial buildings"] <- 1
data$ac[data$treat=="Reference" & # no turbines in industrial landscape
          is.na(data$ac1) & data$ac2=="Industrial buildings"] <- 1

sum(data$ac==1) # Final # of observations
rawdata <- data # save the raw data

data <- data %>% filter(ac==1) # leave attention check passes only

```

## 1.3. Variables

```{r}
## 1) dependent variable
data$support %>% class()

## 2) treatments
data$treat %>% class() # create binary variables indicating the treatment status
data$treat1 <- ifelse(data$treat=="Treatment1", 1, 0); sum(data$treat1==1)
data$treat2 <- ifelse(data$treat=="Treatment2", 1, 0); sum(data$treat2==1)
data$treat3 <- ifelse(data$treat=="Treatment3", 1, 0); sum(data$treat3==1)
data$treat4 <- ifelse(data$treat=="Treatment4", 1, 0); sum(data$treat4==1)
data$treat5 <- ifelse(data$treat=="Treatment5", 1, 0); sum(data$treat5==1)
data$treat6 <- ifelse(data$treat=="Reference", 1, 0); sum(data$treat6==1)
## Seems like there are more attention check failures in the treatments 3-6

## 3) party identification
k <- unique(data$party) %>% length(); unique(data$party) # party4 is the ref
for (i in 1:k) {data[[as.symbol(paste0('party',i))]] <- ifelse(data$party==unique(data$party)[i], 1, 0)}

## 4) gender
k <- unique(data$gender) %>% length(); unique(data$gender) # gender4 is the ref
for (i in 1:k) {data[[as.symbol(paste0('gender',i))]] <- ifelse(data$gender==unique(data$gender)[i], 1, 0)}

## 5) income
k <- unique(data$income) %>% length(); unique(data$income) # income4 is the ref
for (i in 1:k) {data[[as.symbol(paste0('income',i))]] <- ifelse(data$income==unique(data$income)[i], 1, 0)}

## 6) education
k <- unique(data$educ) %>% length(); unique(data$educ) # educ7 is the ref
for (i in 1:k) {data[[as.symbol(paste0('educ',i))]] <- ifelse(data$educ==unique(data$educ)[i], 1, 0)}

## 7) religion
k <- unique(data$relig) %>% length(); unique(data$relig) # relig6 is the ref
for (i in 1:k) {data[[as.symbol(paste0('relig',i))]] <- ifelse(data$relig==unique(data$relig)[i], 1, 0)}

## 8) occupation
k <- unique(data$occu) %>% length(); unique(data$occu) # occu10 is the ref 
for (i in 1:k) {data[[as.symbol(paste0('occu',i))]] <- ifelse(data$occu==unique(data$occu)[i], 1, 0)}

## 9) urban-rural
unique(data$urbanrural)
data$rural[data$urbanrural=="Very rural area"] <- 6
data$rural[data$urbanrural=="Somewhat rural area"] <- 5
data$rural[data$urbanrural=="More rural than urban area"] <- 4
data$rural[data$urbanrural=="More urban than rural area"] <- 3
data$rural[data$urbanrural=="Somewhat urban area"] <- 2
data$rural[data$urbanrural=="Very urban area"] <- 1

## 10) age
data$bornyear <- as.numeric(data$bornyear)
data$age <- 2024-data$bornyear

## 11) reside
data$reside <- as.numeric(data$reside)
```

# 2. Main analysis

## 2.1. Model 1 through 3

This code chunk about is folded for brevity. The reason for using the binary variables not factor ones is to be compatible with the counterfactual analysis in the later sections. 

```{r formulae, class.source="fold-hide"}
# Formula
formula1 <- support ~ treat1 + treat2 + treat3 + treat4 + treat5 + rural +
     party1 + party2 + party3 + age + # party 4 is the ref
     gender1 + gender2 + gender3 + # gender4 is the ref
     income1 + income2 + income3 + income5 + income6 + # income4 is the ref
     educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + # educ7 is the ref
     relig1 + relig2 + relig3 + relig4 + relig5 + relig7 + relig8 + relig9 +
     relig10 + relig11 + # relig6 is the ref
     occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + occu8 + occu9 # occu10 is the ref

formula2 <- support ~ treat1 + treat2 + treat3 + treat4 + treat5 + 
     party1 + party2 + party3 + age + # party 4 is the ref
     gender1 + gender2 + # gender4 is the ref # gender 3 omitted due to singularity
     income1 + income2 + income3 + income5 + income6 + # income4 is the ref
     educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + # educ7 is the ref
     relig1 + relig2 + relig3 + relig4 + relig5 + relig7 + relig8 + relig9 +
     relig10 + relig11 + # relig6 is the ref
     occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + occu8 + occu9 # occu10 is the ref

## formula1_re treats treatment 1 as a reference- for counterfactual scenario purposes
## see 2.2 Section below
formula2_re <- support ~ treat2 + treat3 + treat4 + treat5 + treat6 +
     party1 + party2 + party3 + age + # party 4 is the ref
     gender1 + gender2 + # gender4 is the ref # gender 3 omitted due to singularity
     income1 + income2 + income3 + income5 + income6 + # income4 is the ref
     educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + # educ7 is the ref
     relig1 + relig2 + relig3 + relig4 + relig5 + relig7 + relig8 + relig9 +
     relig10 + relig11 + # relig6 is the ref
     occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + occu8 + occu9 # occu10 is the ref

cfformula2 <- support ~ 1 + treat1 + treat2 + treat3 + treat4 + treat5 + 
     party1 + party2 + party3 + age + # party 4 is the ref
     gender1 + gender2 + # gender4 is the ref # gender 3 omitted due to singularity
     income1 + income2 + income3 + income5 + income6 + # income4 is the ref
     educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + # educ7 is the ref
     relig1 + relig2 + relig3 + relig4 + relig5 + relig7 + relig8 + relig9 +
     relig10 + relig11 + # relig6 is the ref
     occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + occu8 + occu9 # occu10 is the ref

cfformula2_re <- support ~ 1 + treat2 + treat3 + treat4 + treat5 + treat6 +
     party1 + party2 + party3 + age + # party 4 is the ref
     gender1 + gender2 + # gender4 is the ref # gender 3 omitted due to singularity
     income1 + income2 + income3 + income5 + income6 + # income4 is the ref
     educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + # educ7 is the ref
     relig1 + relig2 + relig3 + relig4 + relig5 + relig7 + relig8 + relig9 +
     relig10 + relig11 + # relig6 is the ref
     occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + occu8 + occu9 # occu10 is the ref
```

```{r}
model1 <- lm(formula1, data) # main model
model2 <- lm(formula2, data[data$rural<=3,]) # urban population only, without rural var
model3 <- lm(formula2, data[data$rural>=4,]) # rural population only, without rural var

coeftest(model1, vcov=vcovHC(model1, type="HC1"))[1:10, 1:4] %>% round(digits=2)
coeftest(model2, vcov=vcovHC(model2, type="HC1"))[1:9, 1:4] %>% round(digits=2)
coeftest(model3, vcov=vcovHC(model3, type="HC1"))[1:9, 1:4] %>% round(digits=2)
```

## 2.3. Differences in the expected support across groups and rural-urban

```{r}
# Please manually download two packages below
# https://faculty.washington.edu/cadolph/index.php?page=60
pacman::p_load(simcf, tile)
```

### 1) Urban populations

```{r}
# 1a. comparing sizes, holding constant for farmland (T1 vs T3; T2 vs T3)
pe1a <- coefficients(model2)
vc1a <- vcovHC(model2, type="HC1")
simbetas1a <- MASS::mvrnorm(10000, pe1a, vc1a)

cf_case1a <- cfMake(cfformula2, data[data$rural<=3,], 2)
cf_case1a <- cfChange(cf_case1a, scen=1, "treat1", x=1, xpre=0) # first diff
cf_case1a <- cfChange(cf_case1a, scen=1, "treat2", x=0, xpre=0)
cf_case1a <- cfChange(cf_case1a, scen=1, "treat3", x=0, xpre=1) # first diff
cf_case1a <- cfChange(cf_case1a, scen=1, "treat4", x=0, xpre=0)
cf_case1a <- cfChange(cf_case1a, scen=1, "treat5", x=0, xpre=0)

cf_case1a <- cfChange(cf_case1a, scen=2, "treat1", x=0, xpre=0) 
cf_case1a <- cfChange(cf_case1a, scen=2, "treat2", x=1, xpre=0) # first diff
cf_case1a <- cfChange(cf_case1a, scen=2, "treat3", x=0, xpre=1) # first diff
cf_case1a <- cfChange(cf_case1a, scen=2, "treat4", x=0, xpre=0)
cf_case1a <- cfChange(cf_case1a, scen=2, "treat5", x=0, xpre=0)

cf_result1a <- linearsimfd(cf_case1a, simbetas1a, ci=0.95)

# 1b. comparing sizes, holding constant for industrial (T4 vs T6; T5 vs T6)
model2_re <- lm(formula2_re, data[data$rural<=3,])
pe1b <- coefficients(model2_re)
vc1b <- vcovHC(model2_re, type="HC1")
simbetas1b <- MASS::mvrnorm(10000, pe1b, vc1b)

cf_case1b <- cfMake(cfformula2_re, data[data$rural<=3,], 2)
cf_case1b <- cfChange(cf_case1b, scen=1, "treat2", x=0, xpre=0) 
cf_case1b <- cfChange(cf_case1b, scen=1, "treat3", x=0, xpre=0)
cf_case1b <- cfChange(cf_case1b, scen=1, "treat4", x=1, xpre=0) # first diff
cf_case1b <- cfChange(cf_case1b, scen=1, "treat5", x=0, xpre=0)
cf_case1b <- cfChange(cf_case1b, scen=1, "treat6", x=0, xpre=1) # first diff

cf_case1b <- cfChange(cf_case1b, scen=2, "treat2", x=0, xpre=0) 
cf_case1b <- cfChange(cf_case1b, scen=2, "treat3", x=0, xpre=0) 
cf_case1b <- cfChange(cf_case1b, scen=2, "treat4", x=0, xpre=0) 
cf_case1b <- cfChange(cf_case1b, scen=2, "treat5", x=1, xpre=0) # first diff
cf_case1b <- cfChange(cf_case1b, scen=2, "treat6", x=0, xpre=1) # first diff

cf_result1b <- linearsimfd(cf_case1b, simbetas1b, ci=0.95)

# 1c. comparing landscapes, holding constant for landscapes (T1 vs T4; T2 vs T5)
simbetas1c <- simbetas1a

cf_case1c <- cfMake(cfformula2, data[data$rural<=3,], 2)
cf_case1c <- cfChange(cf_case1c, scen=1, "treat1", x=1, xpre=0) # first diff
cf_case1c <- cfChange(cf_case1c, scen=1, "treat2", x=0, xpre=0)
cf_case1c <- cfChange(cf_case1c, scen=1, "treat3", x=0, xpre=0) 
cf_case1c <- cfChange(cf_case1c, scen=1, "treat4", x=0, xpre=1) # first diff
cf_case1c <- cfChange(cf_case1c, scen=1, "treat5", x=0, xpre=0)

cf_case1c <- cfChange(cf_case1c, scen=2, "treat1", x=0, xpre=0) 
cf_case1c <- cfChange(cf_case1c, scen=2, "treat2", x=1, xpre=0) # first diff
cf_case1c <- cfChange(cf_case1c, scen=2, "treat3", x=0, xpre=0) 
cf_case1c <- cfChange(cf_case1c, scen=2, "treat4", x=0, xpre=0) 
cf_case1c <- cfChange(cf_case1c, scen=2, "treat5", x=0, xpre=1) # first diff

cf_result1c <- linearsimfd(cf_case1c, simbetas1c, ci=0.95)

```


### 2) Rural populations

```{r}
# 2a. comparing sizes, holding constant for farmland (T1 vs T3; T2 vs T3)
pe2a <- coefficients(model3)
vc2a <- vcovHC(model3, type="HC1")
simbetas2a <- MASS::mvrnorm(10000, pe2a, vc2a)

cf_case2a <- cfMake(cfformula2, data[data$rural>=4,], 2)
cf_case2a <- cfChange(cf_case2a, scen=1, "treat1", x=1, xpre=0) # first diff
cf_case2a <- cfChange(cf_case2a, scen=1, "treat2", x=0, xpre=0)
cf_case2a <- cfChange(cf_case2a, scen=1, "treat3", x=0, xpre=1) # first diff
cf_case2a <- cfChange(cf_case2a, scen=1, "treat4", x=0, xpre=0)
cf_case2a <- cfChange(cf_case2a, scen=1, "treat5", x=0, xpre=0)

cf_case2a <- cfChange(cf_case2a, scen=2, "treat1", x=0, xpre=0) 
cf_case2a <- cfChange(cf_case2a, scen=2, "treat2", x=1, xpre=0) # first diff
cf_case2a <- cfChange(cf_case2a, scen=2, "treat3", x=0, xpre=1) # first diff
cf_case2a <- cfChange(cf_case2a, scen=2, "treat4", x=0, xpre=0)
cf_case2a <- cfChange(cf_case2a, scen=2, "treat5", x=0, xpre=0)

cf_result2a <- linearsimfd(cf_case2a, simbetas2a, ci=0.95)

# 2b. comparing sizes, holding constant for industrial (T4 vs T6; T5 vs T6)
model3_re <- lm(formula2_re, data[data$rural>=4,])
pe2b <- coefficients(model3_re)
vc2b <- vcovHC(model3_re, type="HC1")
simbetas2b <- MASS::mvrnorm(10000, pe2b, vc2b)

cf_case2b <- cfMake(cfformula2_re, data[data$rural>=4,], 2)
cf_case2b <- cfChange(cf_case2b, scen=1, "treat2", x=0, xpre=0) 
cf_case2b <- cfChange(cf_case2b, scen=1, "treat3", x=0, xpre=0)
cf_case2b <- cfChange(cf_case2b, scen=1, "treat4", x=1, xpre=0) # first diff
cf_case2b <- cfChange(cf_case2b, scen=1, "treat5", x=0, xpre=0)
cf_case2b <- cfChange(cf_case2b, scen=1, "treat6", x=0, xpre=1) # first diff

cf_case2b <- cfChange(cf_case2b, scen=2, "treat2", x=0, xpre=0) 
cf_case2b <- cfChange(cf_case2b, scen=2, "treat3", x=0, xpre=0) 
cf_case2b <- cfChange(cf_case2b, scen=2, "treat4", x=0, xpre=0) 
cf_case2b <- cfChange(cf_case2b, scen=2, "treat5", x=1, xpre=0) # first diff
cf_case2b <- cfChange(cf_case2b, scen=2, "treat6", x=0, xpre=1) # first diff

cf_result2b <- linearsimfd(cf_case2b, simbetas2b, ci=0.95)

# 2c. comparing landscapes, holding constant for landscapes (T1 vs T4; T2 vs T5)
simbetas2c <- simbetas2a

cf_case2c <- cfMake(cfformula2, data[data$rural>=4,], 2)
cf_case2c <- cfChange(cf_case2c, scen=1, "treat1", x=1, xpre=0) # first diff
cf_case2c <- cfChange(cf_case2c, scen=1, "treat2", x=0, xpre=0)
cf_case2c <- cfChange(cf_case2c, scen=1, "treat3", x=0, xpre=0) 
cf_case2c <- cfChange(cf_case2c, scen=1, "treat4", x=0, xpre=1) # first diff
cf_case2c <- cfChange(cf_case2c, scen=1, "treat5", x=0, xpre=0)

cf_case2c <- cfChange(cf_case2c, scen=2, "treat1", x=0, xpre=0) 
cf_case2c <- cfChange(cf_case2c, scen=2, "treat2", x=1, xpre=0) # first diff
cf_case2c <- cfChange(cf_case2c, scen=2, "treat3", x=0, xpre=0) 
cf_case2c <- cfChange(cf_case2c, scen=2, "treat4", x=0, xpre=0) 
cf_case2c <- cfChange(cf_case2c, scen=2, "treat5", x=0, xpre=1) # first diff

cf_result2c <- linearsimfd(cf_case2c, simbetas2c, ci=0.95)

```

### 3) Plot the expected differences 

```{r}
trace1a <- ropeladder(x=cf_result1a$pe, lower=cf_result1a$lower, upper=cf_result1a$upper,
                      labels=c("Large vs None","Small vs None"),
                      size=0.6, lex=2, lineend="square", plot=1)
trace1b <- ropeladder(x=cf_result1b$pe, lower=cf_result1b$lower, upper=cf_result1b$upper,
                      labels=c("Large vs None","Small vs None"),
                      size=0.6, lex=2, lineend="square", plot=3)
trace1c <- ropeladder(x=cf_result1c$pe, lower=cf_result1c$lower, upper=cf_result1c$upper,
                      labels=c("Large","Small"),
                      size=0.6, lex=2, lineend="square", plot=5)
trace2a <- ropeladder(x=cf_result2a$pe, lower=cf_result2a$lower, upper=cf_result2a$upper,
                      labels=c("Large vs None","Small vs None"),
                      size=0.6, lex=2, lineend="square", plot=2)
trace2b <- ropeladder(x=cf_result2b$pe, lower=cf_result2b$lower, upper=cf_result2b$upper,
                      labels=c("Large vs None","Small vs None"),
                      size=0.6, lex=2, lineend="square", plot=4)
trace2c <- ropeladder(x=cf_result2c$pe, lower=cf_result2c$lower, upper=cf_result2c$upper,
                      labels=c("Large","Small"),
                      size=0.6, lex=2, lineend="square", plot=6)

vertmark1 <- linesTile(x=c(0,0), y=c(0,1), plot=1)
vertmark2 <- linesTile(x=c(0,0), y=c(0,1), plot=2)
vertmark3 <- linesTile(x=c(0,0), y=c(0,1), plot=3)
vertmark4 <- linesTile(x=c(0,0), y=c(0,1), plot=4)
vertmark5 <- linesTile(x=c(0,0), y=c(0,1), plot=5)
vertmark6 <- linesTile(x=c(0,0), y=c(0,1), plot=6)

at.x <- seq(-16, 16, by=4)
tile(trace1a, trace2a, trace1b, trace2b, trace1c, trace2c,
     vertmark1, vertmark2, vertmark3, vertmark4, vertmark5, vertmark6, 
     RxC=c(3,2), limit=c(5,0,0,5),
     xaxis=list(at=at.x),
     height=list(plottitle=3),
     output=list(file="test", width=20/1.618, height=20)
     )

```

## 2.4. Mediation analysis

```{r, warning=F}
## Model 4
data$q7 <- as.numeric(data$q7)
# 7)	Wind turbines are harmful to birds and other wildlife. 
model4_M <- lm(q7 ~ treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
model4_Y <- lm(support ~ q7 + treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)

summary(model4_M)
summary(model4_Y)
model4 <- mediate(model4_M, model4_Y, treat="rural", mediator="q7", boot=T, sims=1000)
summary(model4)

## Model 5
data$q3 <- as.numeric(data$q3)
# 3)	Wind turbines make noises. 
model5_M <- lm(q3 ~ treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
model5_Y <- lm(support ~ q3 + treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
summary(model5_M)
summary(model5_Y)
model5 <- mediate(model5_M, model5_Y, treat="rural", mediator="q3", boot=T, sims=1000)
summary(model5)


## Model 6
data$q6 <- as.numeric(data$q6)
# 6) Wind turbines emit electromagnetic radiation that can cause health-related problems.
model6_M <- lm(q6 ~ treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
model6_Y <- lm(support ~ q6 + treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
summary(model6_M)
summary(model6_Y)
model6 <- mediate(model6_M, model6_Y, treat="rural", mediator="q6", boot=T, sims=1000)
summary(model6)


## Model 7
data$q1 <- as.numeric(data$q1)
# 1)	Addressing climate change should be given high priority, even if it imposes costs on individuals and communities.
model7_M <- lm(q1 ~ treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
model7_Y <- lm(support ~ q1 + treat + rural + 
                 party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
summary(model7_M)
summary(model7_Y)
model7 <- mediate(model7_M, model7_Y, treat="rural", mediator="q1", boot=T, sims=1000)
summary(model7)

## Model 8
data$q2 <- as.numeric(data$q2)
# 2)	Climate change is a more important issue for urban areas than for rural areas.
model8_M <- lm(q2 ~ treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
model8_Y <- lm(support ~ q2 + treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
summary(model8_M)
summary(model8_Y)
model8 <- mediate(model8_M, model8_Y, treat="rural", mediator="q2", boot=T, sims=1000)
summary(model8)

## Model 9
data$q4 <- as.numeric(data$q4)
# 4)	Wind turbines diminish the ruralness of the surrounding area. 
model9_M <- lm(q4 ~ treat + rural + 
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
model9_Y <- lm(support ~ q4 + treat + rural + 
                party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
summary(model9_M)
summary(model9_Y)
model9 <- mediate(model9_M, model9_Y, treat="rural", mediator="q4", boot=T, sims=1000)
summary(model9)

## Model 10
data$q5 <- as.numeric(data$q5)
# 5)	Wind turbines reduce property values. 
model10_M <- lm(q5 ~ treat + rural + 
                 party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
model10_Y <- lm(support ~ q5 + treat + rural +
                  party1 + party2 + party3 +
                  gender1 + gender2 + gender3 + age +
                  income1 + income2 + income3 + income4 + income5  +
                  educ1 + educ2 + educ3 + educ4 + educ5 + educ6 + 
                  relig1 + relig2 + relig3 + relig4 + relig5 + relig6 + relig7 + 
                  relig8 + relig9 + relig10 +
                  occu1 + occu2 + occu3 + occu4 + occu5 + occu6 + occu7 + 
                  occu8 + occu9, data)
summary(model10_M)
summary(model10_Y)
model10 <- mediate(model10_M, model10_Y, treat="rural", mediator="q5", boot=T, sims=1000)
summary(model10)

```

```{r}
coeftest(model4_M, vcov=vcovHC(model4_M, type="HC1"))[7, 1:4] %>% round(digits=2)
coeftest(model4_Y, vcov=vcovHC(model4_Y, type="HC1"))[2, 1:4] %>% round(digits=2)
summary(model4)

coeftest(model5_M, vcov=vcovHC(model5_M, type="HC1"))[7, 1:4] %>% round(digits=2)
coeftest(model5_Y, vcov=vcovHC(model5_Y, type="HC1"))[2, 1:4] %>% round(digits=2)
summary(model5)

coeftest(model6_M, vcov=vcovHC(model6_M, type="HC1"))[7, 1:4] %>% round(digits=2)
coeftest(model6_Y, vcov=vcovHC(model6_Y, type="HC1"))[2, 1:4] %>% round(digits=2)
summary(model6)

coeftest(model7_M, vcov=vcovHC(model7_M, type="HC1"))[7, 1:4] %>% round(digits=2)
coeftest(model7_Y, vcov=vcovHC(model7_Y, type="HC1"))[2, 1:4] %>% round(digits=2)
summary(model7)

coeftest(model8_M, vcov=vcovHC(model8_M, type="HC1"))[7, 1:4] %>% round(digits=2)
coeftest(model8_Y, vcov=vcovHC(model8_Y, type="HC1"))[2, 1:4] %>% round(digits=2)
summary(model8)

coeftest(model9_M, vcov=vcovHC(model9_M, type="HC1"))[7, 1:4] %>% round(digits=2)
coeftest(model9_Y, vcov=vcovHC(model9_Y, type="HC1"))[2, 1:4] %>% round(digits=2)
summary(model9)

coeftest(model10_M, vcov=vcovHC(model10_M, type="HC1"))[7, 1:4] %>% round(digits=2)
coeftest(model10_Y, vcov=vcovHC(model10_Y, type="HC1"))[2, 1:4] %>% round(digits=2)
summary(model10)
```


# 3. Sensitivity analysis

```{r}
sens1 <- medsens(model4, rho.by = 0.01, effect.type = "indirect", sims=1000)
plot(sens1, sens.par="rho", main=" ")
summary(sens1)

sens2 <- medsens(model6, rho.by = 0.01, effect.type = "indirect", sims=1000)
plot(sens2, sens.par="rho", main=" ")
summary(sens2)

sens3 <- medsens(model7, rho.by = 0.01, effect.type = "indirect", sims=1000)
plot(sens3, sens.par="rho", main=" ")
summary(sens3)

sens4 <- medsens(model8, rho.by = 0.01, effect.type = "indirect", sims=1000)
plot(sens4, sens.par="rho", main=" ")
summary(sens4)

sens5 <- medsens(model10, rho.by = 0.01, effect.type = "indirect", sims=1000)
plot(sens5, sens.par="rho", main=" ")
summary(sens5)
```

## Might mediation analysis look different if not using non-parametric?

```{r}
mediate(model4_M, model4_Y, treat="rural", mediator="q7", boot=F) %>% summary()
mediate(model5_M, model5_Y, treat="rural", mediator="q3", boot=F) %>% summary()
mediate(model6_M, model6_Y, treat="rural", mediator="q6", boot=F) %>% summary()
mediate(model7_M, model7_Y, treat="rural", mediator="q1", boot=F) %>% summary()
mediate(model8_M, model8_Y, treat="rural", mediator="q2", boot=F) %>% summary()
mediate(model9_M, model9_Y, treat="rural", mediator="q4", boot=F) %>% summary()
mediate(model10_M, model10_Y, treat="rural", mediator="q5", boot=F) %>% summary()
```

Results are largely the same. 

## PCA analysis

```{r}
data$s1 <- data$q7 # statement 1
data$s2 <- data$q3 # statement 2 (...)
data$s3 <- data$q6
data$s4 <- data$q1
data$s5 <- data$q2
data$s6 <- data$q4
data$s7 <- data$q5

statements <- cbind(data$s1, data$s2, data$s3, data$s4, data$s5,
                    data$s6, data$s7)

facmodel1 <- factanal(statements, factor=3, scores="regression")
print(facmodel1, cutoff=0.0001)

supstate <- cbind(data$support, statements)

facmodel2 <- factanal(supstate, factor=3)
print(facmodel2, cutoff=0.0001)
```

```{r}
round(cor(supstate),2)
```

```{r}
cor.test(facmodel1$scores[,1], data$support)
```



